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Abstract 


This paper explains why high resolution design spaces encourage traditional airfoil opti- 
mization algorithms to generate noisy shape modifications, which lead to inaccurate linear 
predictions of aerodynamic coefficients and potential failure of descent methods. By using 
auxiliary drag constraints for a simultaneous drag reduction at all design points and the 
least shape distortion to achieve the targeted drag reduction, an improved algorithm gen- 
erates relatively smooth optimal airfoils with no severe off-design performance degradation 
over a range of flight conditions, in high resolution design spaces parameterized by cubic 
B-spline functions. Simulation results using FUN2D in Euler flows are included to show 
the capability of the robust aerodynamic shape optimization method over a range of flight 
conditions. 
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chord length of airfoil 
drag coefficient 

gradient of ca with respect to D 
derivative of q with respect to a 
lift coefficient 

gradient of ci with respect to D 

derivative of a with respect to a 

target lift coefficient 

design vector 

mean of random variable 

feasible set for the design vector D 

number of design conditions considered 

free-stream Mach number 

number of design variables 

probability density function of uncertainty parameter u 

number of design conditions 

coordinates of points on plane 

uncertainty parameter 

weighting coefficient 

angle of attack 

average rate of drag reduction at all design points 

minimum rate of drag reduction at all design conditions 

change in design vector 

scalars defining the trust region 

variance of random variable 

a given Mach range 

inner product in Euclidean space 

2-norm in Euclidean space 
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Subscripts and Superscripts 


i index for design condition 

j index for component of design vector 

k index for iteration or iterate 


1 Introduction 

Aerodynamic shape optimization refers to improvement in the performance of an aerody- 
namic surface under geometric shape constraints (such as minimum thickness of an airfoil) 
and performance constraints (such as the target lift requirements). Aerodynamic shape op- 
timization involves computational fluid dynamics (CFD) analysis, sensitivity analysis, and 
nonlinear/nonconvex optimization. Because the performance measures are calculated by 
solving a system of partial differential equations (PDEs), aerodynamic optimization prob- 
lems are PDE-constrained optimization problems. For surveys on aerodynamic shape opti- 
mization and related CFD issues, see references [1]- [4]. In references [1] and [3], Jameson 
examined the use of CFD as a tool for aircraft design. In reference [2], Jameson discussed 
ways to exploit computational simulations more effectively in the overall design process, 
with the primary focus on aerodynamic design, while recognizing the need for an integrated 
multidisciplinary design optimization (MDO). A general overview of industrial applications 
of aerodynamic shape optimization was given by Jameson and Vassberg in reference [4]. 

As computational resources get faster and CFD algorithms are improved, CFD analysis 
and shape optimization will inevitably play a larger role in all phases of the design process. 
These computational procedures will not take the place of wind tunnel and flight-testing 
but can maximize the benefits of expensive experimentation. 

One major obstacle to increased use of high resolution design spaces in aerodynamic 
shape optimization is the existence of many (local) optimal solutions. Jameson et al. 
(ref. [5]) pointed out that “if one considers drag minimization of airfoils in two-dimensional 
inviscid transonic flow, then every shock-free airfoil produces zero drag, and thus optimiza- 
tion based solely on drag has a highly non-unique solution. Different shock-free airfoils can 
be obtained by starting from different initial profiles.” Designers would use the existence 
of multiple solutions to tailor the design to the need of a specific application by using their 
experiences and knowledge of flow physics. On the other hand, automatic numerical algo- 
rithms might be confused by the existence of multiple local optimal solutions. As a result, 
different optimization methods lead to different optimal solutions due to the different search 
paths used in the optimization process. Choices of trust-region and search directions can 
all have noticeable impacts on the final optimal solution. From an aircraft designer’s point 
of view, this dependency of the optimal solution on non-physical numerical tuning param- 
eters is undesirable since the purpose of optimization is to use the underlying flow physics 
to search for a practical design. The dependency of the optimal solution on details of the 
optimization process makes designers question whether the optimization code is guided by 
the physics in the mathematical model or by other unknown factors. One way to avoid these 
questions is to choose a low resolution design space such that the optimal solution in that 
design space is likely unique. In the past, low resolution design space has been constructed 
so that the optimization procedure mimics the manual design process. That is to say, no 
matter what optimization method is used, the optimal solution preferred by designers will 
be recovered. Such a method has limited usefulness since experienced designers can readily 
produce comparable designs without using optimization. 
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Another obstacle to the use of aerodynamic shape optimization is the need to postprocess 
the resulting design. It is known that a high resolution design space tends to generate high 
frequency noise during the airfoil shape modification. Reuther et al. (ref. [6]) mentioned 
that “the use of every surface grid point as an independent design variable also has the 
disadvantage of introducing high-frequency noise into the gradient, that can cause simple 
descent methods to fail.” Therefore, a common practice in airfoil optimization is to filter out 
high frequency components in the search direction during the optimization process (refs. [7]- 
[12]). Both explicit and implicit smoothing procedures for the steepest descent direction were 
used by Szmelter (ref. [13]) for optimization of transonic wings with three design conditions. 
The design principle behind using smoothed search directions is that efficient aerodynamic 
shapes are known to be smooth, as pointed out by Jameson and Vassberg (ref. [4]). In 
practice, an acceptable aerodynamic design must have characteristics that smoothly vary 
with small perturbations in shape and flow conditions. 

To explain why search directions that induce high frequency oscillations must be avoided, 
we will establish a necessary and sufficient condition under which a linear prediction by the 
first-order Taylor expansion is accurate. A mathematical analysis of the gradients for the lift 
and drag shows that a linear prediction (using the linear Taylor expansion) is only accurate 
for smooth variations of the geometry. Such an analysis can explain why optimization 
methods succeed in a low resolution design space. More significantly, the analysis suggests 
a quantitative criterion for an optimization code to check whether a search direction is 
acceptable or not. This analysis leads to various methods for generating reliable descent 
directions even if the design space has very high resolution. 

Another potential problem of using high resolution design space for lift-constrained drag 
minimization is unacceptable off-design performance degradation (ref. [14]). It is well-known 
that lift-constrained drag minimization with a single design point can generate an optimal 
design that is sensitive to variations in the design condition (see, e.g., ref. [15]). A weighted 
average of objective functions corresponding to various design conditions is substituted as the 
objective function to alleviate off-design performance degradation. This substitution leads 
to the standard multipoint optimization model. Unfortunately, the utility of multipoint op- 
timization decreases with an increase in the number of design variables. For example, with 
6 design conditions and 23 design variables (which are the coefficients of sinusoidal basis 
functions), Drela (ref. [14]) discovered that an optimization code tends to exploit insignif- 
icant physical scales during optimization iterations, which leads to off-design performance 
degradation and bumpy optimal airfoil shapes. In other words, given enough resolution in 
the design space, an optimization code might trade an insignificant improvement of aerody- 
namic performance at the design conditions for a severe off-design performance degradation 
(see ref. [16] for a mathematical explanation of how such a phenomenon is related to the 
resolution of the design space) . This behavior is not the fault of the optimization code since 
the off-design performance is not included in the multipoint objective function. 

Off-design performance degradation can be avoided with two approaches: (i) use a well- 
constructed design space so that it is impossible for an optimization code to “over-optimize” 
at the design points, or (ii) treat all design conditions as uncertainties in the design process 
and recast the aerodynamic optimization as an optimization problem under uncertainty 
(refs. [16] and [17]). By treating the design condition as an uncertainty parameter, one 
could explicitly make an optimization code aware of the need to optimize over a range of 
design conditions and help the optimization code to find robust optimal solutions (refs. [16]— 
[18]). 

Finally, an interesting question in optimization research is how to use aerodynamic 
knowledge embedded in flow models to design better aerospace vehicles (that are beyond 
human intuition). To tackle this problem, we must make optimization codes “more intel- 
ligent.” An optimization code can be improved in at least three ways: (i) use high-fidelity 
flow models and/or a more detailed configuration of an aerodynamic vehicle to reduce the 
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error in the knowledge extracted from the mathematical models; (ii) use MDO models to 
make the optimization code aware that the actual aerodynamic performance depends on 
multidisciplinary aspects of the design, such as manufacturing tolerances and wing flexibil- 
ity; and (iii) use optimization under uncertainty to make the optimization code aware that 
the aerodynamic design must adapt to a variety of operating conditions (e.g., high or low 
speeds, with and without flap deflections, fully loaded or empty). Progress in any of these 
research areas will bring MDO for aerospace vehicles closer to reality. 

The main objective of this paper is to show how uncertainty modeling can make an 
optimization code more intelligent in searching for a robust optimal design in a very high 
resolution design space. In section 2, we discuss the importance of smooth airfoil changes 
and the effect of nonsmooth shape modifications on optimization algorithms. In section 
3, we discuss a multiobjective optimization formulation for robust optimization with drag 
minimization as an example. In section 4, we combine the robust optimization strategies 
used in both references [16] and [17] to formulate a random multipoint optimization method 
with auxiliary drag constraints, which provides users with more options for finding a rela- 
tively smooth optimal airfoil with no off-design performance degradation in a high resolution 
design space. The numerical results in section 5 will show how uncertainty modeling, use 
of derivative information, parameterization of the design space, and choice of optimiza- 
tion procedures can all have significant impacts on the final optimal aerodynamic design. 
Concluding remarks are given in section 6. 


2 Accuracy of Linear Prediction of Aerodynamic Coef- 
ficients 

2.1 Accuracy of Sensitivity Derivatives of Aerodynamic Coefficients 

The success of any gradient-based optimization procedure is directly related to the cal- 
culation of accurate sensitivity derivatives. The current state-of-the-art approaches for 
sensitivity analysis for a large number of design variables are based on either discrete or 
continuous adjoint methods. For example, Jameson applied the control theory approach 
to aerodynamic shape optimization that led to adjoint equations for sensitivity calculation 
derived from the potential and Euler equations (ref. [19]). Later, the method was extended 
to Navier-Stokes equations (ref. [20]). The control theory might be applied directly to the 
discrete flow equations that result from the numerical approximation of the flow equations 
by finite element, finite volume, or finite difference methods. This application leads directly 
to a set of discrete adjoint equations with a matrix that is the transpose of the Jacobian ma- 
trix of the full set of discrete nonlinear flow equations (ref. [19]). In theory, if the maximum 
length of the edges in the grid approaches zero, then the difference between the numerical 
derivatives obtained by continuous and discrete adjoint methods also converges to zero (see 
ref. [11] for some numerical confirmation of this hypothesis). Both continuous and discrete 
adjoint methods can generate quite accurate numerical derivatives for aerodynamic shape 
optimization. See references [11] and [12] for a summary on the formulation and develop- 
ment of adjoint methods for sensitivity calculations. Adjoint methods are widely used for 
sensitivity analysis in aerodynamic shape optimization (refs. [6,7,13], and [19]- [38]). 

A sensitivity derivative code is usually verified by comparing its results with the ones 
generated by a finite difference method or a more reliable method based on complex vari- 
ables (see refs. [39] and [40]) for a standard baseline configuration (such as NACA-0012 or 
RAE2822) with a particular parameterization of the boundary curve (such as a B-spline 
curve or a shape function based on Hicks-Henne bump functions). However, we believe that 
a sensitivity derivative code is truly verified when the predicted reduction of the objective 
function based on gradient information is consistently confirmed by the actual reduction of 
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the objective function during the optimization process. Our current and previous numerical 
experiments (refs. [16] and [18]) indicate that the derivative information (with respect to 
shape parameters) at the current airfoil shape given by FUN2D (refs. [21], [41], and [23]) in 
Euler flows can be used reliably to predict the lift and drag coefficients of a slightly modified 
airfoil shape from iteration to iteration. 

To use a gradient-based optimization method to minimize an aerodynamic quantity £, we 
must first understand the conditions under which the derivatives can be used to provide an 
accurate prediction of the value of £ at a new aerodynamic configuration. Mathematically, 
variational analysis is done for shape functions in an appropriate function space (such as a 
Sobolev space). The adjoint equations for derivative calculation derived by Jameson were 
formulated in terms of the variation of the derivative of the boundary transformation that 
maps the boundary of the flow domain to a circle (in two dimensions) or a closed surface 
(in three dimensions) (ref. [19]). Jameson’s formulation of the adjoint equation implies 
that a linear prediction of an aerodynamic quantity is accurate only if the variation of the 
derivative of the boundary transformation is small. 

Suppose that an aerodynamic quantity £ (such as the lift or drag coefficient) is of the 
following form: 

£( r ) = j fr{v) ds, 

where T is the boundary of the domain for the flow solver (such as the airfoil or wing 
surface), fr(v) is defined by the solution of the flow equation, which depends on T, v is the 
vector of coordinates of a point, and the integration (ds) is done with respect to the surface 
area. For simplicity, we only consider the case when T is a two-dimensional curve and has 
the parametric form v = g(t) for 0 < t < 1. Then £(T) can be rewritten as 

£ ( r ) = i) ^ r ( 5 ^) ■ yww dt ’ 

where g'(t) is the derivative of g with respect to t and ||g'(f)|| denotes the length of g'(t). 

Let (r + AT) be a perturbation of the boundary curve with a parametric representation 
v = g(t) + A g(t), 0 < t < 1. Then the variation in £ with respect to AT (or A g) can be 
calculated as follows: 

(Jr A AT) = f 1 [<J r /r,Ar) + {6 v fr,Ag)] • \\g'(t)\\ dt 

+ f fr(g(t))-(g'(t),(Agy(t))-\\g'(t)\r dt, (1) 

Jo 

where (A, B) denotes the value of a linear operator A applied to an element B in an appro- 
priate space. For example, if A and B are vectors in a finite-dimensional Euclidean space, 
then (A,B) denotes the dot product of A and B. In general, we have the following error 
estimate for the first-order sensitivity analysis: 

£(T + AT)-£(T)=(S r £,AT)+o^ (||(A 9 )(t)|| + ||(A fl )'(t)||)dt) . (2) 

From equation (1) we know that the variation of £ with respect to T depends not only 
on the variation AT of T but also on its derivative (A g) 1 . According to equation (2), the 
linear prediction of £ in a neighborhood of T is accurate only if both A g and (A g)' are 
small. Note that (A g)' is small if and only if the length of AT is small. 
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2.2 Importance of Smooth Shape Modification 

The previous analysis sheds new light on the importance of using smooth shape modifications 
during the optimization process. For example, consider a parameterization of the boundary 
curve based on some shape functions: 

m 

v = (po(t) + Yamjt) for 0 < t < 1, 

i= 1 

where represents the baseline shape and tpi(i) are shape functions. By using the 

chain rule, we can get the gradient V 0 £ of £ with respect to a = (ai, . . . , a m ). The linear 
expansion of £ at the point a can be written as 

£(r + AT) - £(T) = <V 0 £, A a) + o(|| Ao||) as ||Aa|| -s- 0. 

However, such an estimate does not provide as much information as equation (2), especially 
when the number of shape functions m is relatively large. 



Figure 1. Curves of sinusoidal function (with m = 22) and its derivative. 

Note that (A g)(t) = Yn=i A a i ■ Vi(t) and (Ag)'(t) = YJi=i A a i ■ is ver y 

large when compared with then a small variation A g does not necessarily mean a 

small (Ag)'. This behavior could happen if the functions ipi are either the sinusoidal basis 
functions or the B-spline basis functions. When m is large, the corresponding sinusoidal 
basis functions are very oscillatory. The oscillation leads to large values of (see fig. 1). 
For the B-spline basis functions, a large value of m means that each basis function is positive 
on a small interval (called the support) and zero everywhere else. This correlation makes a 
B-spline basis function behave like a 5-function with a fast-rising peak near the center of its 
support, which also leads to a large derivative function (see fig. 2). 

Let us use the following numerical example to demonstrate the potential adverse impact 
of the magnitude of (Ag)' on the accuracy of linear predictions of lift and drag. We use 
an airfoil obtained after 50 optimization iterations in reference [18] as the baseline shape 
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Figure 2. Curves of quadratic B-spline function (with 48 knots equally spaced on [0, 1]) and 
its derivative. 
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Figure 3. Accuracy of first-order Taylor expansion with respect to smooth and noisy shape 
modifications. 


because an adverse impact of a large (A g)' can be clearly seen. The airfoil is represented by 
B-splines (as given in fig. 6) with 50 control points. Three control points (one at the leading 
edge and two at the trailing edge) are fixed. The y-coordinates of the remaining 47 control 
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points, denoted by D, and the angle of attack are used as the independent variables. A new 
descent direction AD is calculated with HA-DHoo « 10~ 3 . Then we artificially generate a 
smooth shape modification AD smooth and a noisy shape modification AD noisy as follows: 

(ADsmooth)^ = |(AD)i| and (A£) noisy ) j = (-l)’|(A£)) i | for 1 < i < 47. 

Note that AD smoo th and AD n0 j Sy are identical in terms of the size of their components, but 
ADnoisy represents an oscillatory modification of the airfoil shape and AD smooth represents 
a smooth modification of the airfoil shape. The relative errors of the linear prediction of 
the drag at D are 


e* 



Cd 1 

yD + AD, oti + A a>i, Mij 

1 - 

Cd I 

[■ D,oti,M j) 

- <fg, AS) - Aoifa 

Cd | 

[ D > 

cx$, Mij 

! 


where a* is the angle of attack corresponding to the Mach number Mj. Figure 3(a) shows 
^AD n oi Sy ^ and e, ^AD smoo th) for the linear prediction of drag coefficients. Similarly, figure 

3(b) shows the relative errors of the linear prediction of the lift at D. In both cases, the 
linear prediction of the drag or lift with respect to noisy shape modification is significantly 
less accurate than the linear prediction with respect to smooth shape modification. 


2.3 Conditions for Accurate Linear Prediction 

The previous example can be used to explain why a simple descent method might fail for 
aerodynamic shape optimization with a high resolution design space. Note that if the design 
space has very high resolution, then a numerical approximation of the gradient will contain 
high frequency components. High frequency components have two possible sources: (i) 
numerical approximation errors and (ii) characteristics of the underlying high resolution 
design space. A simple descent method uses a constant multiple of the gradient to modify 
the original airfoil (or wing) shape. Therefore, if AT contains high frequency components 
inherited from the numerical gradient, then (A g)' will be much larger than A g. As a 
consequence, the linear prediction of the value of the objective function at the next iterate 
might be totally wrong (see the residual term in eq. (2)), which could lead to failure of the 
descent method (i.e., the objective function value at the new iterate is actually worse than 
its value at the current iterate). However, if a smoothed gradient is used, then A g does 
not have high frequency noise and (A g)' is as small as A g. Thus, the linear prediction will 
match the actual objective function value at the next iterate, and the descent method will 
operate correctly. 

To demonstrate the importance of using smooth corrections of the airfoil shape during 
the optimization process, we offer the following example. If a noisy airfoil is unintentionally 
generated for some iteration step, it will be very difficult for the optimization code to recover 
from such a mistake. To produce a smooth optimal airfoil and undo the noisy modification 
made in the previous step, A g must be noisy in an equal and opposite direction; however, a 
noisy A g recommended by the linear prediction model might not be reliable. In this case, 
a standard gradient-based optimization method will use the unreliable linear prediction to 
make a shape modification and, as a consequence, either the optimization process fails to 
make any progress or random numerical errors (due to unreliable linear predictions) will 
contaminate the final optimal solution. 

A small (A g)' could also be achieved by using a smoothing procedure for the search 
direction to ensure smooth shape modifications during the optimization process (refs. [7]- 
[13]). Intuitively, larger values of ip\ require a smaller change in coefficients A Gq to keep (A g)' 
at a given level where the linear prediction model is accurate. This requirement gives an 
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incentive for using small step size in gradient-based optimization methods for aerodynamic 
optimization (refs. [11] and [12]). For airfoil parameterization based on B-spline functions, 
we can use a constraint on the variance of Acq , . . . , A a m to achieve smooth modifications 
of airfoils. In fact, for B-spline basis functions, we have J2iLi Witt) = 1 (which is called the 
partition of unit property of the B-spline basis functions (ref. [42])). Thus, i ^(f) = 0. 
As a consequence, a small variance in the data sequence {Aai, . . . , A a m } implies a small 
derivative of A g. Note that one can get a smoothed descent direction by using other data 
smoothing techniques (such as fitting a set of the data by curves with a given magnitude 
of derivatives (ref. [43])). However, we can make the optimization method more intelligent 
by explicitly imposing bound constraints on (A g)' in the trust region formulation for the 
search direction. For example, one could use the following condition to ensure a smooth 
shape modification: 

f ll(Ag)'(f)|| 2 dt < X [ \\Ag(t)\\ 2 dt 
Jo Jo 

or 

max ll(As)'(f)l|oo < A max IKAgX^Hoo, 

where H^Hoo is the supremum norm of a vector v. We can also use the ratio ||(Ag)'||/||Ag|| 
as a quantitative measurement of the quality of the descent direction A g. In general, a 
larger ||(A(?) , ||/|| Ag|| means a less reliable descent direction. 

Later in this paper, we will show that incorporating uncertainty modeling and design 
principles into the optimization procedure can make the procedure smart enough to generate 
relatively smooth descent directions automatically. This approach avoids inaccurate linear 
prediction of aerodynamic coefficients and leads to stable reduction of the objective function 
from iteration to iteration until a local optimal solution is found. 

3 Multiobjective Optimization Formulation for Drag 
Minimization Under Uncertainty 

The main purpose of this section is to provide a multiobjective optimization formulation 
for aerodynamic shape optimization under uncertainty. Mathematically, this formulation 
solves a more difficult problem than the traditional multipoint formulation, but it does not 
require an accurate approximation of the mean and variance of the objective function and 
constraints for its solutions. 

First we would like to point out that optimization under uncertainty is a very broad 
topic, including reliability-based optimization in engineering (ref. [46]), robust optimization 
in engineering design processes (refs. [47]- [49]), and financial optimization (refs. [50] and 
[51]). Recently, optimization under uncertainty was applied in MDO and aerodynamic shape 
optimization (refs. [16]- [18,44,45], and [52]- [55]). 

Some problems in optimization under uncertainty can be formulated as deterministic 
optimization problems. Such formulations seek the best solutions that optimize the objective 
function in the worst case scenario (for the uncertainty parameters) and have been applied 
in truss topology design (refs. [47] and [48]). However, these deterministic formulations are 
known to be NP-hard problems that are computationally intractable (refs. [47] and [48]), 
except the cases when the uncertainty parameters enter as linear terms in the objective 
function or the constraint functions. Since aerodynamic shape optimization problems are 
nonlinear and nonconvex in terms of Mach number M and target lift coefficient c* , such a 
deterministic formulation renders the problem computationally intractable. 

For probabilistic methods to solve optimization problems under uncertainty, the key 
is in the calculation of statistical quantities, such as the mean and variance, and their 
derivatives. The classical robust optimization approach is to minimize both the mean and 
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variance, which is a multiobjective optimization problem with the mean and variance as 
the objective function. In the following we restrict our discussion to a robust optimization 
formulation for aerodynamic shape optimization applications. 

Assume that a vector of uncertainty parameters u is specified for aerodynamic shape op- 
timization. The case of u = M (the Mach number) was discussed in references [16]- [18,44], 
and [45]. Other candidates for uncertainty parameters are Reynold’s number, required lift 
c*, maximum pitching moment and flap setting. Then, the lift-constrained drag minimiza- 
tion under uncertainty can be formulated as the following robust optimization problem: 

zfa") ( £ ( Cd ), C7 ( Q ^) ( 3 ) 


subject to 


D £ T and ci(D,a(u),u) = c^(u) for u £ 0. (4) 

Here c* t (it) is the target lift requirement for the flight condition corresponding to u. T 
is a given feasible set for geometric design variables (that could be defined by geometry 
constraints such as the maximum thickness constraint), and a(u ) is the angle of attack 
corresponding to the uncertainty vector u. The drag and lift coefficients are c d and c/ , 
respectively. The mean and variance of c d are defined as 

E(c d ) = / c d (D,a(u),u) ■ p(u ) du, 

v 2 (cd)= / [c d (D,a(u),u) - E(c d )] 2 p(u) du, 

Jn 

where p(u) is a probability density function of u and fl is the domain of the uncertainty 
vector. 

Unfortunately, equation (3) is only useful to satisfy mathematical curiosity. It is quite 
expensive to compute Ci and c d for any given parameter u. and only crude approximations 
of E(c d ) and a(c d ) can be calculated in practice. It is well-known in the aerodynamic 
optimization community that multiple design conditions can be used in aerodynamic shape 
optimization to produce a shape that is appropriate for a variety of flight conditions. The 
standard approach is to replace fl by a finite sample set and then to find a solution to the 
following multipoint optimization problem: 

r 

min Y^-CdiD, at, ut) (6) 

D,ai,...,a r * — ' 
i= 1 


} ( 5 ) 


subject to 


D G T and ci(D,oti,Ui) = c*^ for 1 < i < r, (7) 

where {ui, . . . , u r } is a set of sample points in f] and Wi are weighting factors that may 
include the effect of p(u). Note that this approach can also be used to solve inverse design 
optimization by using given pressure distributions to guide the optimization method (ref. 
[56]). 

The multipoint optimization approach uses a few flight conditions to deal with off-design 
performance degradation. However, aggressive drag minimization at the design conditions 
could lead to severe off-design performance degradation. In a special case of equation (6) 
where the uncertainty parameter is the Mach number M (i.e., u = M is a scalar), Drela (ref. 
[14]) discovered that if r is much smaller than the number of free design variables, then the 
optimal solution of equation (6) could have severe off-design performance degradation and 
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the optimal airfoil shape could have visible bumps corresponding to the design conditions. 
Drela called this phenomenon point-optimization behavior. A mathematical argument was 
given in reference [16] to explain why it is necessary to use more design points than the 
number of free design variables to avoid undesirable trade-off of marginal improvement at the 
design points and severe off-design performance degradation. However, since computation 
of the aerodynamic coefficients and their derivatives is very expensive, it is not practical to 
use many design points. 

A compromise approach for aerodynamic shape optimization for a range of flight condi- 
tions is to use a few design conditions during the optimization process (so it is computation- 
ally tractable) and to find a Pareto solution of a multiobjective optimization problem while 
trying to avoid off-design performance degradation (refs. [16,18], and [45]). For example, 
to overcome trading marginal improvement at the design points with severe degradation 
of off-design performance when only a few design conditions are used, Huyse and Lewis 
(ref. [17]) and Huyse et al. (ref. [45]) incorporate a random sampling of design points for 
each of the optimization iterations. The motivation for random sampling of design points 
during optimization is to reduce the bias of estimated mean value of the drag due to lack of 
sample size (i.e., the number of design points) and make the optimal solution independent of 
particular set of design points. Li, Huyse, and Padula (ref. [16]) use the profile optimization 
method for drag reduction over the whole range of Mach numbers. The profile optimization 
method solves a sequence of subproblems based on design heuristics. Three design heuristics 
are incorporated in the profile optimization method: (i) use a given linear prediction of drag 
reduction rate to determine the trust region size, (ii) force a simultaneous drag reduction 
at all design points, and (iii) use the least shape distortion to achieve the targeted drag 
reduction. The key idea behind the profile optimization method is to reduce the drag at 
every design point proportionally, which could implicitly reduce the mean and variance of 
the drag over the Mach range. Padula and Li (ref. [18]) provide a unified framework to 
combine these two methods, which leads to a hybrid method including randomly sampled 
design points and the strategies used in the profile optimization method to give the user 
more options of finding a robust airfoil design with respect to the uncertainty in the Mach 
number. 

One unconventional feature of the hybrid method by Padula and Li (ref. [18]) is that 
different objective functions are minimized for each subproblem. This method tries to find 
a locally Pareto optimal solution of the following multiobjective optimization problem: 

min Cd(D, a (u), u) over u £ 0 (8) 

D,a(u) 

subject to the constraints given in equation (4). 

Recall that (. D*,a*(u )) is locally Pareto optimal for (8) if ( D*,a*(u )) satisfies (4) and 
there is no ( D,a{u )) in a neighborhood of ( D*,a*(u ) such that ( D,a(u )) satisfies (4) and 
Cd(D,a(u),u) < Cd{D*,a*(u),u) for all u £ with Cd(D,a(u),u) < Cd(D* ,a* {u),u) for at 
least one u in 0. If D* is a strict local optimal solution to (6), let a*(u) be the angle of attack 
such that ci(D* , a* (u) , u) = c*(u) for u € fi, then (D* , a* (u)) is also locally Pareto optimal 
for (8). Note that by changing the objective function during optimization iterations, we are 
not solving a particular single objective optimization problem. Instead, such an approach is 
designed to find a local Pareto optimal solution of (8). See the next section for more details. 

4 The Randomized Multipoint Optimization Method 
With Auxiliary Drag Constraints 

In this section, we introduce the constrained multipoint optimization method for finding a 
Pareto optimal solution of equation (8) . To facilitate our presentation, we divide this section 
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into two subsections: the first one contains the description of our main algorithm and the 
second one provides a detailed explanation of the technical background for each step in the 
algorithm. 


4.1 Algorithm Description 

Our main algorithm is based on a randomized multipoint optimization formulation with 
auxiliary drag constraints. This algorithm uses a random sampling of design points and 
smart descent directions to find an approximate solution of equation ( 3 ). 

The Randomized Multipoint Optimization Method (RMOM) With Auxiliary 
Drag Constraints. Let D° be a given initial design vector, e > 0 a tolerance control for 
termination, 77 > 0 random sampling control parameters, 7 m i n the minimum reduction rate 
for each design point, 7 & ve(> 7min) the reduction rate for the weighted average of the drag 
at all design points, and k = 0 . Then construct a sequence of design vectors as follows: 

1 . Generate random samples of design points. Find r design vectors ui jk , ■ ■ ■ , u r 7 (when 
k > 2) in the uncertainty region such that 


■^i,k 


^i,k — 1 


< 77 for 1 < i < r. 


2. Compute feasible angles of attack. Find 077, 0:27, ■ ■ ■ , such that 


ci{D k ,ati, k ,Ui, k ) = c* Li for 1 < i < r. 

3 . Perform adaptive adjustment of the integration scheme. Let Y^l= 1 W 'i-M ' c d{D, on, Ui tk ) 
be a numerical approximation of the expected drag E(cd ) given in equation ( 5 ), where 
wi, k ,W2, k , • ■ • , w r , k are positive weights that include the effect of the probability den- 
sity function p(u) . 

4 . Check the termination criterion. Consider the following two conditions: 


^ ) 10 i, k 

i = 1 


< 9 c d , n *. , 

-Q^{D h ,a itk ,u itk ) 


< e 


( 9 ) 


and 


where 


mm 




i= 1 


^ Tmin 

< C 

OWe 


dc d , nt: , 

-Q^{D h ,a itk ,u itk ) = 


dc d ( 
dD ( 


da, J OCi 

dD 


(D , Q77 , . 


(10) 


If either equation ( 9 ) or equation ( 10 ) is satisfied, then output D k as an optimal design 
vector and terminate the algorithm. 


5 . Solve a trust region subproblem. Let c d ,i, k and cu.k be the linear approximations of 
the drag and lift at ( D k ,oti jk ): 


/ o \ o 

ci,i, k (&D,Aai) = ci(D k ,a itk ,u itk ) + 1 -^,Ad\ + 
Cd,i,k(AD, Aon) = c d (D k ,a itk ,u itk ) + 1 -^,Ad\ + 
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where the derivatives are evaluated at (D k , ck*,* , it*,*). Consider the following trust 
region subproblem: 


min V' w itk ■ c dt i :k (AD, Aa t ) subject to (11) 

AD,Acxi ' 
i= 1 

— If i,kfa < A Di < j it kSk for 1 < i < m, 

-a itk < A oti < a max - a itk for 1 < i < r, 

a }i}k (AD, A oti) = for 1 < i < r, 

Cd,i,kA D A*i) < (1 -7min) ■ Cd(D k ,a it k ,u itk ) for 1 < i < r, 

where 7 i tk > 0 are some scaling factors. Choose 5 k > 0 such that equation ( 11 ) is 
feasible and the least norm solution (A D k , Actu^, . . . , Aa r>k ) of equation ( 11 ) satisfies 
the following condition: 

r r 

^ ' Cd,i,k(AD ' , Ao^fc) — (1 7ave) ^ ) '^i,k ' ^‘d f ^ ^i, k ) • 

i= 1 i= 1 

6 . Generate the new iterate. Let oti jk+ 1 = + A ai tk for 1 < i < r and D k+1 = 

D k + A D k . 

7. Start a new iteration. Update k by k + 1 and go back to step 1. 

4.2 Technical Discussion of Algorithm 

The previously described algorithm has a few unconventional strategies adopted from ref- 
erences [16] and [45] for finding robust optimal solutions. The algorithm is based on robust 
design heuristics for aerodynamic optimization. To fully understand how the algorithm finds 
a robust optimal solution, we need to explain what happens at each step in the algorithm. 




Figure 4. Two optimal airfoils generated by the RMOM with fixed and random design 
points. 

Step 1 implements random sampling of the design points from iteration to iteration. 
This sampling allows an optimization code to modify the airfoil shape at various locations 
on the airfoil that influence performance to various design conditions and, thus, eliminate 
the point-optimization behavior. To illustrate, we plot the optimal airfoils generated by 
the RMOM with fixed design points and random design points, respectively, in figure 4. 
Note that the random sampling of the design points makes the optimal solution smoother. 
Figure 4(b) shows a close-up view of the boxed area in figure 4(a). This example shows 
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some advantage of having the option of using random sampling of the design points in step 
1 to make the optimal solution less dependent on the choice of design points. The random 
sampling control parameter t*, can change from iteration to iteration or be set to zero if 
fixed design points are preferred. 

In step 2 we find the feasible angles of attack. When the target lift is less than the 
maximum achievable lift, the lift coefficient is always a monotone increasing function of 
the angle of attack in a neighborhood of Therefore, the optimal choice of the angle 
of attack is always the smallest one that satisfies the lift constraint. Thus, the inequality 
constraint for the lift in equation (3) can be replaced by an equality constraint to generate 
the same optimal solution. Kim et al. (ref. [34]) modified a standard flow solver to adjust 
the angle of attack during the iteration process so that the flow solver could calculate the 
angle of attack for a given target lift with a computational time only marginally more than 
the calculation of the lift and drag for a given angle of attack. 

In step 3, one can choose arbitrary weights to define the objective function. In terms of 
multiobjective optimization for the drag at the given design points, the weights determine 
a particular weak optimal solution on the weak efficient front of the Pareto optimization 
problem. However, the choice of weights based on the probability density function p(u) and 
a numerical integration formula for the expected drag makes the approach related to the 
minimization of the expected value of the drag. 

In step 4, the termination criterion given is of theoretical nature since it is impossible 
in practice to determine what should be an appropriate value of e. However, this criterion 
can serve as an assessment tool to measure the quality of the iterates generated by the 
algorithm. 

The two termination conditions can be derived by using the implicit differentiation rule. 
Here we are assuming that, for any flight condition and a feasible target lift, one can adjust 
the angle of attack to achieve the target lift. For each design vector D, we can find the 
smallest oti(D) such that 

ci(D,ai(D),Ui) = c*j. (12) 

Assuming that T is the whole space, we can rewrite equation (6) as 

mm 1(D) (13) 


where 

r 

!(D ) = £>,- c d (D,ai(D),Ui). 
i = 1 

Using the implicit differentiation, we obtain 


dl 

dD 




( 

1 dD 



(14) 


Note that equation (14) was discussed by Reuther et al. (ref. [6]) using variational analysis 
approach and later incorporated by Kim et al. (ref. [34]) into optimization algorithms for 
solving multipoint optimization problems. Thus, the first-order optimality condition for the 
multipoint optimization problem is 



This condition explains why we use equation (9) as a stopping criterion. 
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min 


= 0 . 


(16) 


The first-order optimality condition for a minimax version of the multipoint optimization 
problem is 

f, ( Be a (fg) 

' l dD (d £L '\ 

1=1 \ \dai ) 

Note that if equation (16) is satisfied, then it is impossible to find a descent direction that 
will reduce the drag at all design points simultaneously. 

If there is no requirement on a minimum reduction of the drag at each design point 
(i.e. , 7 m i n < 0), then equation (10) never holds and the RMOM is essentially searching 
for a multipoint optimal solution with respect to dynamically allocated design points. If 
7min = 7ave, then equation (9) implies equation (10) and the RMOM will be terminated by 
equation (10). In this case, the RMOM generates an optimal solution in a minimax sense. 
In general, the ratio provides an option for choosing a combination of the weighted 
average and the minimax strategies. 

In step 5, the linear subproblem formulated is similar to the linear subproblem for a 
trust region method based on sequential linear programming, except for the auxiliary drag 
constraints. The auxiliary constraints ensure a simultaneous predicted drag reduction at 
all design points, when 7 min > 0. If 7 m i n is a large negative number (say, —1000) and the 
design points are fixed from iteration to iteration, then the RMOM solves the multipoint 
optimization problem because the drag constraints in equation (11) become superfluous 
and are always satisfied. Similarly, if 7 m i n is negative but the design points are randomly 
sampled from iteration to iteration, then the algorithm becomes a variation of the expected 
value optimization method (refs. [17] and [45]). If we choose 7 m i n = 7 ave and fix the design 
points from iteration to iteration, then the RMOM is reduced to the profile optimization 
method proposed by Li, Huyse, and Padula (ref. [16]). Therefore, by using different values 
of 7min and sampling schemes in the RMOM, we can get various hybrids of the multipoint 
optimization method, the expected value optimization method, and the profile optimization 
method. 

In step 6, the new iterate generated is chosen to increase reliability of the linear prediction 
and to avoid random shape distortions. For example, let S be the orthogonal complement 
of the linear space generated by the gradients of the lift and the drag with respect to the 
design variables; that is, if 


S = <v 


j Qq \ 

(v, g^(D k ,ai :k ,u itk )\ = 0 and 

(v, |^( D k ,a itk ,u itk )^ = 0 for 1 < i < r| 


then vectors in S do not change the linear prediction of the lift and drag at the design 

, £ 

points. Thus, if AH is in 5, then for any real number t, (A D k + tAD ) is also a solution 

of the trust region subproblem as long as (A D k + tAD ) is in the trust region. In the case 
when the number of design variables is greater than 2 r, S is the set of solutions of 2 r linear 
equations with more than 2 r variables. Thus, the dimension of S is at least one and the 
linear subproblem might have multiple solutions. The simplex method finds a vertex of the 
solution region, and the interior point method finds the analytical center of the solution 
region. Therefore, different solvers of the linear subproblem might generate different values 
of A D k , which lead to different iterates. These differences explain the impact a solver of 
the linear subproblem could have on the intermediate iterates. If a quasi-Newton method 
is used to find A D k , then the solution of the trust region subproblem is unique. In this 
case, the S'-components in the search direction A D k are not introduced randomly. Instead, 
they reflect the (random) numerical error of the quadratic approximation of the objective 
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function. Note that any shape change in S is unnecessary and could be considered as a 
random perturbation in the descent direction, which might lead to random shape distortions 
during the optimization iterations. From previous studies (e.g., refs. [7]- [13]), we know that 
smoothing of the airfoil shape is a quite common procedure added in an airfoil optimization 
process to filter out high frequency components in the airfoil shape to prevent noises in 
airfoils. To avoid randomly adding elements in 5 to the descent direction AD, we force AD 
to be the least norm solution of the trust region subproblem, which has two advantages: 
(i) it generates a solution AD k that is orthogonal to vectors in S, which prevents random 
additions of those vectors in AD k , and (ii) the descent direction AD k has the smallest norm 
among all solutions of the trust region subproblem, which helps to increase the reliability 
of the linear prediction of the lift and drag. 

The new iterate generated in step 6 is influenced by the values of drag reduction rate 
parameters 7 m i n and 7 ave . When 7 m j n > 0, the purpose of the auxiliary drag constraints is 
to find a common descent direction for the drag at the design points. To demonstrate the 
implication of such a strategy, we eliminate ai in equation (11) and rewrite the auxiliary 
drag constraints as follows: 

AD , (D , OL^ k , ^ 7minCrf(-fl , OL^ k , U^fc) for 1 < i < T. 

Multiplying both sides of the above inequality by > 0 and adding all the resulting 
inequalities together, we obtain 

(\AD k ,Y^^i^(D k ,a itk ,u itk )\j < - ^i'yminCd(D k ,a itk ,u itk )^j <0 

for \i > 0 and i ^ i > 0- If the gradient of the drag at an off-design point u is a 
nonnegative combination of the gradients of the drag at the design points, then AD k is also 
a descent direction for the drag at u. For example, one can prove that A D k is a descent 
direction for the drag at u whenever is in a neighborhood V of the convex cone generated 
by the gradients of the drag at the design points. In general, a larger 7 m i n implies a larger 
neighborhood V ; that is, increasing 7 m i n reduces the risk of trading a drag reduction at the 
design points with off-design drag increment. 

In the special case when 7 m i n = 7 ave , our algorithm is actually solving a discrete version 
of the classical robust optimization problem given by equation (3). In fact, let E(cd) (or 
E{cd)) and a 2 (cd) (or a 2 (cd)) denote the discrete version of the expected value and variance 
of Cd at D k + AD k (or D k ) based on the sample points Ui jk ; that is, 

r 

B{C(i) — ^ ^ ^i,k ’ 5 A 

i= 1 

r 2 

° 2 ( C d) = (c d ,i,k( A OLi, k ) - E{c d )) , 

i= 1 
r 

E( c d) = ^2 w hk -c d (D k ,ai, k ,Ui, k ), 
i= 1 

r 

d 2 (cd) = ( Cd(D k ,ati, k ,Ui , k ) - E(c d )) . 

i= 1 

In most circumstances, if the design space has enough resolution, then one could find a 
descent direction such that the equalities hold for the auxiliary drag constraints in equation 
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(11). In this case, we have 

E(c d ) = (1 - 7 min )E(c d ) and a 2 (c d ) = (1 - 7 min) 2 h 2 (c d ). (17) 

That is, the predicted mean and variance of the drag at the design points are reduced 
by the given rate 7 m in- If the linear prediction is accurate, then the algorithm actually 
minimizes both estimates of the mean and variance of the drag computed at the design 
points. Therefore, the algorithm generates a solution of a discrete version of the classical 
robust optimization problem given by equation (3). 

Finally, note that the change of each design variable is proportional to the scaling factors 
7 i t k- One common choice for 7 ^ is the absolute value \D%\ of the ith component of D k . 
However, other choices of 7^ might be necessary for particular applications. Scaling factors 
7 i t k have a significant impact on how the optimization code will work. In general, a given 
percentage of predicted drag reduction can be achieved in many ways. By using various 
values of 7 , one can force the optimization code to focus on modifying certain parts of 
the aerodynamic shape to achieve the predicted drag reduction. 

Theoretically, we could use the standard trust region strategy to check the ratio of the 
actual drag reduction and the predicted drag reduction and update 7 m i n and 7 ave accord- 
ingly. We conjecture that the behavior of such an algorithm is expected to be similar to 
traditional trust region algorithms: an actual sufficient reduction of the objective function 
from iteration to iteration and any cluster point of the iterates being a stationary point of 
either a multipoint optimization problem or a minimax optimization problem. 

Note that any local strict Pareto optimal solution of the following multiobjective opti- 
mization problem, 

min c d (D,ai,Ui ) over 1 < i < m, 

is also a Pareto optimal solution of equation (8). Thus, a Pareto optimal solution of equation 
(8) can be found by using a few design points. In practice, we usually choose j ave to achieve 
a meaningful drag reduction (such as 1-4 percent) . The choice of 7 m i n is based on robustness 
requirement as mentioned previously. To avoid technical complications, we do not include 
the adjustment of j ave and 7 m i n in the algorithm, but readers could experiment with different 
values of j ave and 7 m i n for different iterations. 

In summary, the RMOM is aimed at finding a Pareto optimal solution, while minimizing 
an approximation of the expected drag and maintaining a level of simultaneous drag reduc- 
tion at all design points from iteration to iteration. The descent direction is selected in such 
a way that a given predicted drag reduction rate is achieved with a minimum amount of 
shape modification. This strategy helps to eliminate high frequency components in the de- 
scent direction, and reduces the need for airfoil smoothing during the optimization iterations 
(even if the airfoil is represented by B-splines with many free design variables) . 

5 Numerical Simulation Results 

All our numerical results indicate that the RMOM generates relatively smooth airfoils with 
no severe off-design performance degradation, even when the design space consists of B- 
spline curves with over 50 control points. In this section, we only include a few sample 
results to illustrate how the sampling of design points, resolution of the design space, and 
the choice of minimum drag reduction rate (i.e., the value of 7 m i n ) affect the optimal solution 
generated by the RMOM. 

5.1 Technical Information for Simulation Setting 

For a given airfoil, the lift and drag coefficients q and c d are obtained by solving Euler 
equations using FUN2D (refs. [21], [41], and [23]). For numerical results presented in this 
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Figure 5. Grid used to solve Euler equations by FUN2D. 



Figure 6. Parameterization of NACA-0012 by cubic B-splines with 50 control points. 

paper, an unstructured grid with 2872 points and 5744 elements is used (see fig. 5 for the 
grid near the airfoil). The derivatives of the lift and drag are obtained by solving the adjoint 
equations of the discrete flow equations (i.e., by using the discrete adjoint method). With 
the given grid, the flow solver and the adjoint solver are terminated when the 2-norms of 
the residual of the density equations and its adjoint counterpart are reduced by at least 
5 orders of magnitude. The initial airfoil is NACA-0012. The airfoil is parameterized by 
using B-splines with either 23 or 50 control points. See figure 6 for a parameterization of the 
NACA-0012 airfoil with 50 control points. The x-coordinates of all the control points are 
fixed during optimization. However, the maximum thickness of the airfoil is not constrained. 
The y-coordinates of the first and last control points are fixed such that the trailing edge of 
the airfoil is fixed during optimization. The y-coordinate of the leading edge of the airfoil is 
also fixed such that the chord length of the airfoil remains constant during optimization. The 
number of free design variables is the total number of control points minus 3. The weights 
7 i t k in the RMOM are chosen to be \D k | (i.e., proportional to the size of design variables). 
For simplicity, we consider only the lift-constrained drag minimization with uncertain flight 
speed (i.e., u = M) and q = c* = 0.4 for all i. For each iterate D k , we find angles of 
attack Qi^, . . . , a TOj fc such that \ci{D k , a^k, Ui t k) — c* J/c* 4 < 0.001. (Note that Kim et al. 
(ref. [34]) replaced the equality lift constraint c; = c] by 0 < ( c* — c*)/c* < 0.003.) The 
probability density function p{u) is the uniform distribution over the Mach range [0.7, 0.8]. 
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The results presented in this section are selected from the following three sets of test cases: 

1. Four equally spaced design points in [0.7, 0.8], 47 free design variables, 7 ave = 3 percent, 
and 7 m i n = —100000, 1, 2, or 3 percent 

2. Four randomly distributed design points in [0.7, 0.8], 47 free design variables, 7 ave = 3 
percent, and 7 m j n = —100000, 1, 2, or 3 percent 

3. Four equally spaced design points in [0.7, 0.8], 20 free design variables, j ave = 3 percent, 
and 7min = —100000 or 3 percent 

All the optimization runs are terminated after either 50 or 100 iterations, which allows us 
to compare different test cases easily. In this paper, we did not use the derived stopping 
criteria to terminate the optimization runs. 

5.2 Optimal Airfoil Shapes 

The optimal airfoils generated by the RMOM are quite smooth and do not have four bumps 
corresponding to four design conditions (i.e., Mach numbers). This is different from the 
numerical simulation results obtained earlier by Drela on multipoint airfoil optimization 
when sinusoidal basis functions were used to parameterized the airfoil shape. Figure 7 is a 
typical optimal airfoil generated by the RMOM after 50 iterations. All the optimal airfoils 
have shape patterns that are surprisingly similar to the airfoils generated by a two-point 
optimization method (Mi = 0.8, c* Ll = 0.3438 and M 2 = 0.82 , c /,2 = 0.3348), starting from 
NACA-0012 (see figs. 2 and 6 in ref. [26]). Interestingly, the RMOM actually finds such a 
shape pattern to be optimal among all possible airfoil shapes represented by B-splines with 
47 free design variables. In comparison, the airfoil generated by Elliott and Peraire (ref. [26]) 
was given in the form of z(x) = zo(x ) + J2i=i Pi^i( x ), where hi(x) are specific Hicks-Henne 
bump functions that allow the airfoil to deform to a highly aft-loaded supercritical airfoil. 
Thus, we contend that the RMOM has a potential of finding relatively realistic airfoil designs 
without intervention from aerodynamic experts. 



Figure 7. Typical airfoil generated by the RMOM after 50 iterations starting from baseline 
airfoil NACA-0012. 


5.3 Impact of Sampling of Design Points 

Our numerical results show that random sampling of design points can be an effective means 
to prevent the optimizer from excessively modifying any particular location on the airfoil 


19 




surface to reduce the shock corresponding to a particular flight condition. This sampling 
avoids specific shock bumps on the airfoil surface corresponding to design conditions, which 
could occur during the standard multipoint airfoil shape optimization (ref. [14]). As a 
consequence, one might get smoother airfoil shapes by using randomized design points (see 
fig. 4). In terms of drag reduction, whether we use random design points or fixed design 
points, the drag profiles after 50 iterations are quite similar (see fig. 8). Figure 9 shows the 
complete history of sampled design points for 50 iterations. One can see that the variations 
of design points from iteration to iteration are quite substantial. 



Figure 8. Drag curves versus Mach numbers for 50th iterates generated by the RMOM with 
7min = —100000 percent, 7 ave = 3 percent, and 47 free design variables. 



Iteration Number 


Figure 9. Complete history of sampled design points for 50 iterations of the RMOM with 
7min = —100000 percent, 7 ave = 3 percent, and 47 free design variables. 

From a theoretical point of view, unless the size of random perturbation of design points 
(i.e. , 7 *,) is sufficiently small, we can not ensure that the new iterate is indeed an “im- 
provement” of the old iterate because changing design points means changing the objective 
function. Practically, changing design points could increase computational expenses. In 
general, for fixed design points, we do not have to adjust the angles of attack to satisfy the 
lift constraints because the linear predictions of the lift are accurate enough for the new 
iterate to satisfy the lift constraints. In contrast, if randomized design points are used, then 
the angles of attack must be adjusted to satisfy the lift constraints in each iteration. For the 
three sets of test cases reported here, each iteration requires 4 function calls and 8 derivative 
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calls to FUN2D for the RMOM with fixed design points, while the RMOM with randomized 
design points requires approximately 8 additional function calls to find the feasible angles 
of attack at the 4 new design points. However, if the flow solver has an option of producing 
the drag and the required angle of attack for a given lift input, then the computational cost 
of using randomized design points is the same as using fixed design points. 

5.4 Impact of Auxiliary Drag Constraints 

Auxiliary constraints only become meaningful if 7 m i n > 0, which allows one to enforce a 
certain percent of predicted drag reduction at each design point. By varying 7 m i n from — oo 
to 7ave, we hope to achieve a dynamic balance between the most conservative optimization 
strategy (with 7 m i n = 7 ave ) and the most aggressive optimization strategy (with 7 m i n = 
— oo). Figure 10 includes the drag curves for the solutions generated by the RMOM after 
50 iterations using various values of 7 m i n . All these iterates were generated with 7 ave = 3 
percent and 47 free design variables. In general, the drag curve for the solution has a smaller 
slope if 7mi n is larger. Using the mean and standard deviation of the drag at 24 equally 
spaced sample Mach numbers in the Mach range [0.7, 0.8] as approximations of the actual 
mean and standard deviation of the drag on [0.7, 0.8], we attained the results presented in 
table 1. Table 1 clearly shows that 7 m i n provides an option for trade-off between the mean 
and standard deviation of the drag for the optimal solution over the given Mach range. 



Figure 10. Drag curves versus Mach numbers for solutions generated by the RMOM after 
50 iterations using various values of 7 m in- 


Table 1. Estimated Mean and Standard Deviation of Drag in Mach Range [0.7, 0.8] 



Estimated Mean 

Estimated Standard Deviation 

Omin — OO 

0.0038 

0.0025 

7min = 2 percent 

0.0039 

0.0023 

7min = 3 percent 

0.0040 

0.0018 


5.5 Impact of Resolution of Design Space 

Figure 11 shows the the drag curves versus Mach numbers for the 100th iterate generated by 
the RMOM with 7 m i n = 3 percent and 7 ave = 3 percent for 20 and 47 free design variables, 
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respectively. For Mach numbers near 0.8, the drag for the 100th iterate using 20 free design 
variables is twice as large as the drag for the 100th iterate using 47 free design variables. 
By checking the optimization history, we find that the drag coefficients start to oscillate 
after 50 iterations when using 20 free design variables. Figure 12 shows the drag coefficients 
for 4 fixed design points from the 51st to the 100th iteration. In other words, if the initial 
baseline is the 50th iterate, then the RMOM could not improve the initial design with the 
most conservative optimization strategy (using 7 m i n = 7 ave ) because the design space does 
not have enough resolution to provide a highly restrictive search direction for drag reduction 
at each of the design points. However, if a different optimization strategy is employed (with 
7min = — 100000 percent < < 7 ave = 3 percent) , then it is much easier to find a less restrictive 
search direction for an average drag reduction at the given design points. As a consequence, 
the drag curve of the 100th iterate generated by using 20 free design variables is not much 
different from that of the 100th iterate generated by using 47 free design variables (see fig. 
13). 



Figure 11. Drag curves versus Mach numbers for 100th iterates generated by the RMOM 
with 7 min = 3 percent and 7 ave = 3 percent. 


x 10 3 



Iteration Number 


Figure 12. Complete history of drag values at 4 fixed design points for iterations 51-100. 
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Figure 13. Drag curves versus Mach numbers for 100th iterates generated by the RMOM 
with 7 min = —100000 percent and j ave = 3 percent. 

6 Concluding Remarks 

This research is motivated by the needs of conceptual and preliminary design teams. These 
teams desire smooth airfoil shapes that are similar to the baseline design but have improved 
drag performance over a range of flight conditions. The resulting airfoil optimization method 
modifies a large number of design variables to search for the true optimal solution, yet it 
avoids any severe off-design performance degradation that is characteristic of single-point 
methods. 

New versions of commercial aircraft are variants of baseline designs with modified op- 
erating conditions such as required lift or cruise speed. Given a good initial design, this 
new optimization method improves the performance without making drastic shape changes. 
Small reductions in the expected value of drag translate into large reductions in fuel con- 
sumption. Moreover, the design team gains valuable information by examining the details 
of the flow field around the airfoil before and after optimization. 

As a demonstration of the optimization algorithm, we studied a lift-constrained drag 
minimization problem for a two-dimensional airfoil in Euler flow with 47 B-spline coefficients 
as free design variables. An unstructured grid computational fluid dynamics code, FUN2D, 
predicts the lift and drag values and their gradients with respect to changes in airfoil shape 
and changes in angle-of-attack. We were able to generate relatively smooth airfoils with no 
severe off-design performance degradation over the range of Mach numbers, even when the 
number of design points was as small as four. 

This work further demonstrates that it is possible to avoid off-design performance degra- 
dation by using robust airfoil shape optimization methods, even if a large number of free-form 
design variables are used. Post-optimization analysis shows that the optimized airfoil also 
has better performance than the original airfoil at 24 off-design conditions. Results shown in 
the paper are typical of all numerical experiments performed by the authors. It is important 
to note that the solution of each optimization subproblem is a relatively smooth airfoil and 
that each airfoil has improved performance over the previous one. Termination criteria are 
provided and these criteria help the user judge when to stop because further improvement 
is unlikely. 

To understand the importance of using smooth shape modification for aerodynamic 
shape optimization, we investigate the linear Taylor approximation error of an integral of 
the pressure around the airfoil. Such an error study shows that assuming the exact gradient 
information is given, the accuracy of the linear Taylor approximation of the lift or drag 
coefficient depends not only on the magnitude of shape change, but also on the magnitude 
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of the derivative of shape change. As a consequence, with the same order of magnitude 
of shape modification, noisier shape modifications (i.e., larger derivatives of shape change) 
lead to larger disparities between the predicted reduction and the actual reduction of the 
drag. Thus, noisy shape modifications for a descent method may increase the drag for 
lift-constrained drag minimization problems. 

Even though we do not use airfoil smoothing during the optimization, the least norm 
solution of the trust region subproblem is implemented to reduce random airfoil shape 
distortions during the optimization process, which forces a minimal shape change to achieve 
a predetermined rate of reduction for the given objective function. 

Note that randomized multipoint optimization method (RMOM) does not solve an op- 
timization problem with a single objective function. Using different objective functions 
for formulating optimization subproblems allows us to enforce a robust optimization policy 
without any information on the mean and variance of the drag. RMOM attempts to find a 
local Pareto optimal solution of the multiobjective optimization problem (8) that behaves 
like an optimal solution to the robust optimization problem (3). 

The algorithm presented in our paper allows us to choose a combination of the weighted 
average and minimax strategies, which corresponds to a trade-off of performance and robust- 
ness during the optimization process, to find better aerodynamic shapes in a high resolution 
design space. 
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